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One-dimensional detrended fluctuation analysis (ID DFA) and multifractal detrended fluctuation 
analysis (ID MF-DFA) are widely used in the scaling analysis of fractal and multifractal time series 
because of being accurate and easy to implement. In this paper we generalize the one-dimensional 
DFA and MF-DFA to higher-dimensional versions. The generalization works well when tested 
with synthetic surfaces including fractional Brownian surfaces and multifractal surfaces. The two- 
dimensional MF-DFA is also adopted to analyze two images from nature and experiment and nice 
scaling laws are unraveled. 
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I. INTRODUCTION 

Fractals and multifractals are ubiquitous in natural 
and social sciences 0. The most usual records of ob- 
servable quantities are in the form of time series and their 
fractal and multifractal properties have been extensively 
investigated. There are many methods proposed for this 
purpose Syl, such as spectral analysis, rescaled range 
analysis 0,11 H B Ej , fluctuation analysis [T(i| , de- 

], wavelet trans- 



trended fluctuation analysis [Til ^ 
form module maxima (WTMM) E 



trended moving average 0, El l2lll22tl23lj~to list a few 
It is now the common consensus that DFA and WTMM 
have the highest precision in the scaling analysis Q>H[24|. 

The idea of DFA was invented originally to investi- 
gate the long-range dependence in coding and noncoding 
DNA nucleotides sequence Then it was generalized 
to study the multifractal nature hidden in time series, 
termed as multifractal DFA (MF-DFA) [d|. Due to the 
simplicity in implementation, the DFA is now becoming 
the most important method in the field. 

Although the WTMM method seems a little bit com- 
plicated, it is no doubt a very powerful method, es- 
pecially for high-dimensional objects, such as images, 
scalar and vector fields of three-dimensional turbulence 
EE El Elm El. In contrast, the original DFA method 
is not designed for such purpose. In a recent paper, a 
first effort is taken to apply DFA to study the roughness 
features of texture images [3(|. Specifically, the DFA is 
applied to extract Hurst indices of the one-dimensional 
sequences at different image orientations and their aver- 
age scaling exponent is estimated. Unfortunately, this is 
nevertheless a one-dimensional DFA method. 

In this work, we generalize the DFA (and MF-DFA as 
well) method from one-dimensional to high-dimensional. 
The generalized methods are tested by synthetic surfaces 
(fractional Brownian surfaces and multifractal surfaces) 
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with known fractal and multifractal properties. The nu- 
merical results are in excellent agreement with the the- 
oretical properties. We then apply these methods to 
practical examples. We argue that there are tremendous 
potential applications of the generalized DFA to many 
objects, such as the roughness of fracture surfaces, land- 
scapes, clouds, three-dimensional temperature fields and 
concentration fields, turbulence velocity fields. 

The paper is organized as follows. In Sec. [HI we rep- 
resent the algorithm of the two-dimensional detrended 
fluctuation analysis and the two-dimensional multifrac- 
tal detrended fluctuation analysis. Section ITTll shows the 
results of the numerical simulations, which are compared 
with theoretical properties. Applications to practical ex- 
amples are illustrated in Sec. IIVI We discuss and con- 
clude in Sec. 



II. METHOD 

A. Two-dimensional DFA 

Being a direct generalization, the higher-dimensional 
DFA and MF-DFA have quite similar procedures as the 
one-dimensional DFA. We shall focus on two-dimensional 
space and the generalization to higher-dimensional is 
straightforward. The two-dimensional DFA consists of 
the following steps. 

Step 1: Consider a self-similar (or sclf-affinc) surface, 
which is denoted by a two-dimensional array X(i,j), 
where i = 1, 2, • • • , M, and j = 1, 2, • • • , N. The sur- 
face is partitioned into M s x N s disjoint square seg- 
ments of the same size s x s, where M s = [M/s] and 
N s = [N/s]. Each segment can be denoted by X VyW such 
that X v ^ w (i,j) — X(li + j) for 1 ^ i,j ^ s, where 

li = (v — l)s and 1% = (w — l)s. 

Step 2: For each segment X v w identified by v and w, 
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the cumulative sum u v>w (i,j) is calculated as follows: 



B. Two-dimensional MF-DFA 



fe,=l k 2 =i 



(1) 



where 1 ^ i, j ^ s. Note that u v . w itself is a surface. 

Step 3: The trend of the constructed surface u VfW can 
be determined by fitting it with a prcchosen bivariate 
polynomial function u. The simplest function could be a 
plane. In this work, we shall adopt the following five de- 
trending functions to test the validation of the methods: 



^V,W (^3 j ) 
^"U,W (^3 j) 
Uy^yj (i, j) 
Uy^yj (i, j) 



ai + + C 

„-2 , 1-2 , 



(2) 
(3) 
(4) 
(5) 

ai 2 + bj 2 + cij + di + ej + / , (6) 



at' + bj~ + c 
aij + bi + cj + d , 
ai 2 + bj 2 + ci + dj + e 



where 1 ^ i,j ^ s and a, 6, c, d, e, and / are free pa- 
rameters to be determined. These parameters can be es- 
timated easily through simple matrix operations, derived 
from the least squares method. We can then obtain the 
residual matrix: 



(7) 



The detrended fluctuation function F (v, w, s) of the seg- 
ment X vw is defined via the sample variance of the resid- 
ual matrix e v , w (i,j) as follows: 



s s 

F 2 (v,w,s) = — ^2^2^v,w(iJ) 



(8) 



i=l o = l 



Note that the mean of the residual is zero due to the 
dctrcnding procedure. 

Step 4: The overall detrended fluctuation is calculated 
by averaging over all the segments, that is, 



F 2 (s) 



1 



M s N s 



M,JV« 



^^F 2 ( V , W , S ). (9) 



V— 1 w= 1 



Step 5: Varying the value of s in the range from s m ; n w 
6 to s max ~ min(Af, N)/4, we can determine the scaling 
relation between the detrended fluctuation function F(s) 
and the size scale s, which reads 



F(s) ~ s H 



(10) 



where H is the Hurst index of the surface @, El HE j 
which can be related to the fractal dimension by D = 

3-h HQ. 

Since N and M need not be a multiple of the segment 
size s, two orthogonal trips at the end of the profile may 
remain. In order to take these ending parts of the surface 
into consideration, the same partitioning procedure can 
be repeated starting from the other three corners |3lT | . 



Analogous to the generalization of one-dimensional 
DFA to one-dimensional MF-DFA, the two-dimensional 
MF-DFA can be ascribed similarly, such that the 
two-dimensional DFA serves as a special case of the 
two-dimensional MF-DFA. The two-dimensional MF- 
DFA follows the same first three steps as in the two- 
dimensional DFA and has two revised steps. 

Divide a self-similar (or self-affine) surface X(i,j) into 
M s x N s (M s = [M/s] and N s = [N/s]) disjoint phalanx 
segments. In each segment X v>w (i,j) compute the cumu- 
lative sum u(i,j,s) using Eq. QJ. With one of the five 
regression equations, we can obtain u(i,j, s) to represent 
the trend in each segment, then we obtain the fluctuate 
function F(v,w,s) by Eq. (JSJ). 

Step 4: The overall detrended fluctuation is calculated 
by averaging over all the segments, that is, 



F q (s) = 



1 



1/9 



M„N, 



(11) 



1 ui— 1 



where q can take any real value except for q = 0. When 
q = 0, we have 

r 1 Ms Ns } 

^ v= 1 u>— 1 ) 

according to L'Hospital's rule. 

Step 5: Varying the value of s in the range from s m j n « 
6 to s max w min(M, N)/4, we can determine the scaling 
relation between the detrended fluctuation function F q (s) 
and the size scale s, which reads 



F q (s) ~ s 



h{q) 



(13) 



For each q, we can get the corresponding traditional 
r(q) function through 



r(q) = qh(q) - D f 



(14) 



where Df is the fractal dimension of the geometric sup- 
port of the multifractal measure H i s thus easy to 
obtain the generalized dimensions D q j^M HE IH3 and 
the singularity strength function a(q) and the multifrac- 
tal spectrum /(a) via Legendre transform |38| |. In this 
work, the numerical and real multifractals have Df = 2. 
For fractional Brownian surfaces with a Hurst index H, 
we have h(q) = H, 



C. A note on the generalization 

To the best of our knowledge, the first few steps of the 
one-dimensional DFA and MF-DFA in literature are or- 
ganized in the following order: Construct the cumulative 
sum of the time series and then partition it into segments 
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of the same scale without overlapping. In this way, a di- 
rect generalization to higher dimensional space should be 
the following: 

Step I: Construct the cumulative sum 

i 3 

fc 1= i k 2 =i 

Step II: Partition u(i,j) into N s x M s disjoint square 
segments. The ensuing steps are the same as those de- 
scribed in Sec. Ill Al and Sec. Ill Bl 

It is easy to show that, for the one-dimensional DFA 
and MF-DFA, the residual matrix in a given segment 
is the same no matter which step is processed first, 
either the cumulative summation or the partitioning. 
This means that we have two manners of generalizing to 
higher-dimensional space, that is, Steps 1-2 in Sec. Ill Al 
and Steps I-II aforementioned. Our numerical simula- 
tions show that both these two kinds of generalization 
gives the correct Hurst index for fractional Brownian 
surfaces when adopting two-dimensional DFA. However, 
the two-dimensional MF-DFA with Steps I-II gives wrong 
r(q) function for two-dimensional multifractals with an- 
alytic solutions where the power-law scaling is absent, 
while the generalization with Steps 1-2 does a nice job. 

The difference of the two generalization methods be- 
comes clear when we compare u VlW (i,j) in Eq. Q with 
u(i,j) in Eq. lfT3|) . We see that u V)W (h +i, l<z +j) is local- 
ized to the segment X v , w , while u(l\ + i, 1% + j) contains 
extra information outside the segment when i < l\ and 
j < I2, which is not constant for different i and j and 
thus can not be removed by the detrending procedure. 
In the following sections, we shall therefore concentrate 
on the correct generalization expressed in Sec. Ill Al and 
Sec.lHBl 

III. NUMERICAL SIMULATIONS 
A. Synthetic fractional Brownian surfaces 

We test the two-dimensional DFA with synthetic frac- 
tional Brownian surfaces. There are many different meth- 
ods to create fractal surfaces, based on Fourier transform 
filtering [34 . l39j |. midpoint displacement and its vari- 
ants 111 to 1 4JJ , ci rculant embedding of covariance ma- 
trix [43, |43 L 4J, |45j , periodic embedding and fast Fourier 
transform 46j| . top-down hierarchical model 01 > an d so 
on. In this paper, we use the free MATLAB software Fr- 
acLab 2.03 developed by INRIA to synthesize fractional 
Brownian surfaces with Hurst index H . 

In our test, we have investigated fractional Brown- 
ian surfaces with different Hurst indices H ranging from 
0.05 to 0.95 with an increment of 0.05. The size of the 
simulated surfaces is 500 x 500. For each H, we gen- 
erated 500 surfaces. Each surface is analyzed by the 
two-dimensional DFA with the five bivariate functions 
in Eqs. IL'Kil) . The results are shown in Fig. ^ We can 



see that the estimated Hurst indexes H are very close 
to the preset values in general. The deviation of Hurst 
index H becomes larger for large values of H . 




H 



FIG. 1: Comparison of the estimated Hurst index H using 
Eqs. 12161 with the true value H. The error bars show the 
standard deviation of the 500 estimated H values. The re- 
sults corresponding to Eqs. il-'illil are translated vertically for 
clarity. 



In Fig. |2 we show the log-log plot of the detrended 
fluctuation F (s) as a function of s for two synthetic frac- 
tional Brownian surfaces with H = 0.2 and H = 0.8, 
respectively. There is no doubt that the power-law scal- 
ing between F(s) and s is very evident and sound. Hence, 
the two-dimensional DFA is able to well capture the self- 
similar nature of the fractional Brownian surfaces and 
results in precise estimation of the Hurst index. 

We also adopted fractional Brownian surfaces to test 
the two-dimensional multifractal detrended fluctuation 
analysis. Specifically, we have simulated three frac- 
tional Brownian surfaces with Hurst indexes Hi = 0.2, 
H2 = 0.5, and H3 = 0.8, respectively. The five regression 
equations are used in the detrending. We calculated 
h(q) for q ranging from —10 to 10 according to Eq. (T3J. 
All the F q (s) functions exhibit excellent power-law scal- 
ing with respect to the scale s. The function r(g) can 
be determined according to Eq. i|14f> . The resultant r(q) 
functions are plotted in Fig. |21 with the inset showing the 
h(q) functions. We can find from the figure that, for each 
surface, the five functions of r(q) (and h(q) as well) cor- 
responding to the five detrending functions collapse on a 
single curve. Moreover, it is evident that h(q) = H and 
r(q) = qH — 2. The three analytic straight lines intersect 
at the same point (q = 0, r(g) = —2). These results are 
expected according to theoretical analysis. 

We stress that, when fractional Brownian surfaces are 
under investigation, both the two-dimensional DFA and 
MF-DFA can produce the same correct results even when 
Steps I-II arc adopted. 
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FIG. 2: Log-log plots of the detrended fluctuation function 
F(s) with respect to the scale s for H — 0.2 (top panel) and 
H = 0.8 (bottom panel) using Eqs. (J2EU- The lines are the 
least squares fits to the data. The results corresponding to 
Eqs. are translated vertically for clarity. 



B. Synthetic two-dimensional multifractals 

Now we turn to test the MF-DFA method with syn- 
thetic two-dimensional multifractal measures. There ex- 
its several methods for the synthesis of two-dimensional 
multifractal measures or multifractal rough surfaces |26j| . 
The most classic method follows a multiplicative cascad- 
ing process, which can be either deterministic or stochas- 
tic |4J,|49|,|5fi|5l| The simplest one is the p-model pro- 
posed to mimick the kinetic energy dissipation field in 
fully developed turbulence [49l |. Starting from a square, 
one partitions it into four sub-squares of the same size 
and chooses randomly two of them to assign the measure 
of p/2 and the remaining two of (1 — p)/2. This parti- 
tioning and redistribution process repeats and we obtain 
a singular measure p. A straightforward derivation fol- 
lowing the partition function method |3Sj results in the 
analytic expression: 



T{q)=q-l-\0g 2 [p1 + (l-p) q ] 



(16) 



A relevant method is the fractionally integrated sin- 
gular cascade (FISC) method, which was proposed to 
model multifractal geophysical fields [52j and turbulent 
fields The FISC method consists of a straightfor- 




FIG. 3: Plots of r(q) extracted by using the five detrending 
functions 12I6H as a function of q. The three straight lines are 
r(q) = qH - 2 for Hi = 0.2, Hi = 0.5, and H 3 = 0.8, respec- 
tively. The inset shows the corresponding h(q) functions. 



ward filtering in Fourier space via fractional integration 
of a singular multifractal measure generated with some 
multiplicative cascade process so that the multifractal 
measure is transformed into a smoother multifractal sur- 
face: 



f(x) = fi(x) <g> \x\ 



-(1-H) 



(17) 



where g) is the convolution operator and H £ (0, 1) is 
the order of the fractional integration [2(|, whose r(q) 
function is IMIHl: 



(18) 



r(q) = q(l + H) - 1 - log 2 [p« + (1 - p) q ] 



The third one is called the random W cascade method 
which generates multifractal rough surfaces from random 
cascade process on separable wavelet orthogonal basis 

St- 

In our test, we adopted the first method for the syn- 
thesis of two-dimensional multifractal measure. Starting 
from a square, one partitions it into four sub-squares of 
the same size and assigns four given proportions of mea- 
sure pi = 0.05, P2 = 0.15, P3 = 0.20, and P4 = 0.60 to 
them. Then each sub-square is divided into four smaller 
squares and the measure is redistributed in the same way. 
This procedure is repeated 10 times and we generate mul- 
tifractal "surfaces" of size 1024x1024. The resultant r(q) 
functions estimated from the two-dimensional MF-DFA 
method are plotted in Fig. where the inset shows the 
h(q) functions. We can find that the five functions of r(q) 
(and h(q) as well) corresponding to the five detrending 
functions collapse on a single curve, which is in excellent 
agreement with the theoretical formula: 



r(q) = - log 2 (pi + p q 2 + p\ + pi) 



(19) 



We stress that, when we use Steps I-II instead of Steps 
1-2, the resultant r(q) estimated by the MF-DFA method 
deviates remarkably from the theoretical formula. In- 
deed, the power-law scaling for most q values is absent 
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FIG. 4: Plots of r(q) extracted by using the five detrending 
functions 12I6H as a function of q. The continuous line is the 
theoretical formula l|19|l . The inset shows the corresponding 
h(q) functions. 



and thus the alternative algorithm with Steps I-II and the 
resulting r(q) is completely wrong. In addition, we see 
that different detrending functions give almost the same 
results. The linear function @ is preferred in practice, 
since it requires the least computational time among the 
five. 



IV. EXAMPLES OF IMAGE ANALYSIS 
A. The data 

In this section we apply the generalized method to an- 
alyze two real images, as shown in Fig. [5] Both pictures 
are investigated by the MF-DFA approach since it con- 
tains automatically the DFA analysis. The first example 
is the landscape image of the Mars Yardangs region [3(| , 
which can be found at |http: / / sse.jpl.nasa.gov| The size 
of the landscale image is 2048 x 1536 pixels. The second 
example is a typical scanning electron microscope (SEM) 
picture of the surface of a polyurethane sample foamed 
with supercritical carbon dioxide. The size of the foam- 
ing surface picture is 1200 x 800 pixels. 




FIG. 5: Left: The image of the Yardangs region on the Mars. 
Right: A scanning electron microscope picture of the surface 
of a polyurethane sample foamed with supercritical carbon 
dioxide. 



The SEM picture of the surface of a polyurethane sam- 
ple were prepared in an experiment of polymer foaming 
with supercritical carbon dioxide. At the beginning of the 
experiment, several prepared polyurethane samples were 
placed in a high-pressure vessel full of supercritical car- 
bon dioxide at saturation temperature for gas sorption. 
After the samples were saturated with supercritical CO2, 
the carbon dioxide was quickly released from the high- 
pressure vessel. Then the foamed polyurethane samples 
were put into cool water to stabilize the structure cells. 
Pictures of the foamed samples were taken by a scanning 
electron microscope. 

The two images were stored in the computer as two- 
dimensional arrays in 256 prey levels. We used Eq. |j2J 
for the detrending procedure. The two-dimensional ar- 
rays were investigated by the multifractal detrended fluc- 
tuation analysis. For each picture, we obtained the r(g) 
function and the h(q) function as well. If r(q) is nonlin- 
ear with respect to q or, in other words, h(q) is dependent 
of q, then the investigated picture has the nature of mul- 
tifractality. 



B. Analyzing the Mars landscape image 

We first analyze the Mars landscape image shown in 
the left panel of Fig.EJwith MF-DFA. FigureEJilhistrates 
the dependence of the detrended fluctuation F q (s) as a 
function of the scale s for different values of q marked 
with different symbols. The continuous curves are the 
best linear fits. The perfect collapse of the data points 
on the linear lines indicates the evident power law scal- 
ing between F q (s) and s, which means that the Mars 
landscape is self-similar. 
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FIG. 6: Log-log plots of the detrended fluctuation function 
F q (s) versus the lag scale s for five different values of q. The 
continuous lines are the best fits to the data. The plots for 
q = —3, q = 0, q = 3, and q = 6 are shifted upwards for 
clarity. 



The slopes of the straight lines in Fig. give the es- 
timates of h(q) and the function r(g) can be calculated 
accordingly. In Fig. is shown the dependence of r(q) 
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with respect to q for — 6 ^ q ^ 6. We observe that r(q) is 
linear with respect to q. This excellent linearity of r(q) is 
consistent with the fact that h(q) is almost independent 
of q, as shown in the inset. Hence, the Mars landscape 
image does not possess multifractal nature. 
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FIG. 7: Dependence of r(q) with respect to q. The solid line 
is the least squares fit to the data. The inset plots h(q) as a 
function of q. 



The values of h(q) are estimated by the slopes of the 
straight lines illustrated in Fig.|Hlfor different values of q. 
The corresponding function r(g) is determined according 
to Eq. I|14|) . In Fig. [5] is illustrated r(q) as a function of 
q for — 6 ^ q ^ 6. We observe that r(q) is nonlinear with 
respect to q, which is further confirmed by the fact that 
h(q) is dependent of q, as shown in the inset. The non- 
linearity of r(g) and h(q) shows that the foaming surface 
has multifractal nature. 























2.4 


©« 












S 2 .2 
2 


0. 

* 




J 






















1.8 


-5 


5 
9 






-4 


-2 





2 


4 


< 



9 

FIG. 9: Dependence of r(q) with respect to q. The inset 
shows h(q) as a function of q. 



C. Analyzing the foaming surface image 



Similarly, we analyzed the foaming surface shown in 
the right panel of Fig. with the MF-DFA method. Fig- 
ure |HI illustrates the dependence of the detrended fluc- 
tuation F q (s) as a function of the scale s for different 
values of q marked with different symbols. The continu- 
ous curves are the best linear fits. The perfect collapse 
of the data points on the linear line indicates the evident 
power law scaling between F q (s) and s, which means that 
the Foaming surface is self-similar. 
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FIG. 8: Loglog plots of the detrended fluctuation function 
F q (s) versus the lag scale s for five different values of q. The 
continuous lines are the best fits to the data. The plots for 
q — —3, q — 0, q = 3, and q = 6 are shifted upwards for 
clarity. 



V. DISCUSSION AND CONCLUSION 

In summary, we have generalized the one-dimensional 
detrended fluctuation analysis and multifractal de- 
trended fluctuation analysis to two-dimensional versions. 
Further generalization to higher dimensions is straight- 
forward. We have found that the higher-dimensional 
DFA methods should be performed locally in the sense 
that the cumulative summation should be conducted af- 
ter the partitioning of the higher-dimensional multifrac- 
tal object. Extensive numerical simulations validate our 
generalization. The two-dimensional MF-DFA is applied 
to the analysis of the Mars landscape image and foam- 
ing surface image. The Mars landscape is found to be a 
fractal while the foaming surface exhibits multifractality. 

At last, we would like to stress that there are tremen- 
dous potential applications of the generalized DFA in 
the analysis of fractals and multifractals. In the two- 
dimensional case, the methods can be adopted to the 
investigation of the roughness of fracture surfaces, land- 
scapes, clouds, and many other images possessing self- 
similar properties. In the case of three dimensions, it 
could be utilized to qualify the multifractal nature of 
temperature fields and concentration fields. Possible ex- 
amples in higher dimensions are stranger attractors in 
nonlinear dynamics. Concrete applications will be re- 
ported elsewhere in the future publications. 
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